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Abstract 


The injectivity tests consist of injecting a phase, usually water, through wells into an oil reservoir in 
order to collect information about its parameters. Knowing these reservoir’s parameters can be valuable in 
order to improve oil production. In that sense, many studies concerning injectivity tests in multilayered 
reservoirs have been presented. However, these studies regard only commingled systems, that is, layers 
that are separated by impermeable barriers. Therefore, the present work developed an analytical model in 
the Laplace space for multilayered reservoirs under injectivity tests considering formation crossflow along 
all its layers. The accuracy of the proposed solution was verified by comparison with a finite difference 
flow simulator. The results provided by the analytical model and by the numerical data were consistently 
similar. Furthermore, the data obtained by the analytical solution was used to estimate the reservoir’s 
equivalent permeability. Calculated values presented a satisfactory accuracy for all cases. 


1. Introduction 


An injectivity test consists on injecting water into an oil reservoir and, in response to that, many 
significant information about the reservoir can be inferred, which can be valuable in order to improve the oil 
production. For instance, equivalent permeability, outer boundary condition and recoverable oil volume can 
be inferred. The way of detecting those pieces of information is through analyzing the pressure transient 
response. In this work, an analytical solution is proposed concerning multilayered reservoirs with formation 
crossflow. 


Formation crossflow, different than wellbore crossflow, happens when there is a nonzero vertical perme- 
ability at the interface between two adjacent layers. In this work formation crossflow was defined like in the 
work of Sun and Gao (2017), which presented a semipermeable-wall system. Besides formation crossflow, 
this work will consider that all layers properties, such as porosity, permeability and thickness, may be 
different in each layer. 


Based on Sun and Gao (2017), the work proposed by Ehlig-Economides and Joseph (1987), presented 
an analytical model for a multilayered reservoir under single phase flow considering formation crossflow. 
Properties may also be different from layer to layer. In addition, the work presented by Ehlig-Economides 
and Joseph (1987), was developed in the Laplace field and used numerical inversion of Laplace transforms, 
more specifically, the Stehfest Algorithm (Stehfest, 1970) in order to find a pressure response in the real 
field. The present work will also be presented in the Laplace field and the Stehfest algorithm will be used 
in order to convert it into the real field. 


The formulation proposed in this work, combines the formation crossflow feature with a two phase 
flow model. To model the two phase flow composition, this work was based on the one presented by 
LF Bittencourt Neto et al. (2020b) which extended the solution from a model for a radially composed 
reservoir under single phase flow as in given by LF Bittencourt Neto et al. (2020a), because during injectivity 
tests in a homogeneous reservoir, there are still different regions, but in respect to the fluids. The formulation 
proposed by LF Bittencourt Neto et al. (2020a) was also done in the Laplace domain, considered different 
properties in a two-layered model. Furthermore, in LF Bittencourt Neto et al. (2020b), a single layered 
reservoir was considered and formation crossflow was unconsidered. 

The work proposed by Barreto et al. (2011), presented a model for injectivity tests like in the work of 
Bratvold et al. (1990), however Barreto et al. (2011) considered a multilayered reservoir. The solution 
proposed by Barreto et al. (2011) for the injection period assumes, like the present work, a radially infinite 
reservoir. Both the works proposed by Bratvold et al. (1990) and Barreto et al. (2011) did not considered 
formation crossflow. 

This present work is innovating as it presents a formulation considering formation crossflow for multilayered 
reservoirs under injectivity tests. This paper will be divided, essentially, in description of the model, 
mathematical formulation, results and discussions and conclusions. 
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2. Model Description 


Before the injectivity test begins, the reservoir is filled with single phase solution (oil), then, the injection 
of water in the oil reservoir begins, establishing the two phase fluid flow scenario. For this case, consider an 
infinite, multilayered, homogeneous and isotropic reservoir. In addition, it has an oil viscosity (j1,) that 
is constant, as well as constant water viscosity (Ww). It also has constant and small total compressibility 
(cz). The porosity (¢;) may have a different value in each layer, as well as the thickness (h;). Constant 
flow rate is considered. Formation damage and wellbore storage will be disregarded. Moreover, consider a 
piston-like displacement, i.e., the water front in each layer advances uniformly. All the formulation was 
done in consistent units. 

Even though the rock properties within each layer are the same, there are different composition properties 
within them, such as viscosity. Hence, this two phase flow model can be represented by a two regions model 
with moving boundary, with the water region increasing with time. Because of their similarities, the single 
phase model for composite reservoirs can be extended to model the two phase case LF Bittencourt Neto 
et al. (2020b). 

Constant layer flow rates and thicknesses will be considered and there is formation crossflow within these 
layers. In each layer, region 1 will represent the region that contains the water being injected and region 2 
the region that contains oil only. Figure 1 represents that system: 





Fig. 1 — Two phase fluid flow in a reservoir 


In figure 1, r,, is the radius of the injector well, ry the water front radius, which will be defined later 
for each time t and r, is the external radius. Consider Ap,;; = p;; — p;, the pressure change in each layer 


j =1,2,...,n is governed by the following diffusion equation for each region of fluid i = w,o: 
2 OAD ji 
(kh) jV"Apji = (bh) jeri + (Agi — Apg—ayr) Xj-1 — (Apgar — Anji) X; [1] 


The coefficient of semipermeability X,; is given as in Ehlig-Economides and Joseph (1987) for layers 7 
and j + 1 relating the permeabilities of each one of this layers: 
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2 
~ Q(Ahg) /Rug] + x54 + x3 





Xj [2] 


where X, = Xp = 0. Ah; and k,, are the thickness and the permeability of shale between layers j and 
j +1 and x; is defined as: 


h 2 
ie [3] 
a5 


where k,, is the vertical permeability of layer 7. If there is no formation crossflow between layers j and 
j +1 then, X; is zero. 


Observing the reservoir laterally, figure 2 illustrates this model for a two layered case: 


hy Layer 1 


hg Layer 2 





ray 122 


Fig. 2 — Two phase fluid flow in a two layer reservoir with formation crossflow 


Notice that, since a single region of permeability is considered per layer, there is only a single semiperme- 
ability coefficient in each adjacency. Eq.(2) evidences that only the rock properties are considered which 
are the same for all layers. If there were more regions in respect to the permeability, there would be more 
semipermeability coefficients. 

Now, consider the volume injected in layer 7 (Barreto et al., 2011): 


v,(t) = / qj(t)dt (4) 
0 


To calculate the water front radius (rf) for each layer j the formula given by Buckley et al. (1942) for 
cylindrical geometry is considered: 
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By 
rip(t) = (rw? + sengehy tus 5] 


where Tig represents the derivative of the curve of fractional flow of water at layer j at each point of 
fluid saturation. Since a piston flow is considered, f/,; is given by (LF Bittencourt Neto et al., 2020b): 


1 


i oe a 
Soe 


Here S¥,, is the oil residual saturation and Si. is the irreducible water saturation in layer 7. The water 
and oil fluid mobility of layer 7 (\7), which will be used in the initial boundary condition, are defined as: 


sj kirk 


nN , where 7 = w,o (7] 
Li 





a 


Here ki is the permeability of layer 7 related to phase in region 7 and kj; is layer 7 permeability. 


3. Mathematical Formulation 


The formulation presented in this section was mainly extended from the work by Ehlig-Economides and 
Joseph (1987) for the formation crossflow consideration and from the work by LF Bittencourt Neto et al. 
(2020a) for the extension of a radial composite reservoir for a two phase formulation. 


Considering the hypotheses from the previous section, the model description is given by the diffusion 
equation [1] and its suitable boundary conditions: 


e The Initial Condition (I.C.) which occurs at time t = 0 and reflects that the reservoir is initially in 
balance, that is, the pressure is the same in all layers (apart from the hydrostatic effect). 


e The Internal Boundary Condition (I.B.C) which is related to the fluid injection in each layer. 


e The External Boundary Condition (E.B.C) which refers to the flow behavior at the extreme limit of 
the reservoir and, in this work, an infinite-like reservoir is considered. 
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Consider Kj; = (kjhj) and w; = (¢;h;) and let i index each region of phase F' = w, 0, then the following 
equations model the problems for all layers 7 = 1,...,n and regions i = 1,2 (Lefkovits et al., 1961; 
Ehlig-Economides and Joseph, 1987): 


Region 1 (rw <r < rjF): 








OAp, (r,t 
PDE: «jV?Apji = Lefty PED + (Apji — Apg—1y1) Xj-1 — (Apgaai — Apji)X53,t >0 [8] 
IC: Apyi(r,0) =0 [9] 
_ “ OAp;i 
IBC: qi= —2niwh(r ) - [10] 
Region 2 (r > rjr): 
OAp,ja(r,t 
PDE : KjyV? Apya = wje4ptg PED + (Apjo — Apg—1)2)Xj-1 — (Apgaiye — Apja)Xj3,t>0 [11] 
IC: Apje(r,0) =0 [12] 
EBC: lim Apjo(r = 00, t) = 0 [13] 


When j = 1, X;-1 = 0 and when j = n, X; = 0. Consider p;, to be the Laplace transform of the pressure 
variation and u to be the Laplace’s variable, then, in the Laplace field, equations [8] to [13] become: 


Region 1: 


_ 6 OAD 51 





[15] 


T=Tw 


Region 2: 


ODE : Kj VAD jo = Wi CLUMoAD 59 (7, t) + (ADj2 — ADG—1)2)Xj—-1 — (APG 412 — ADj2) Xj; 7 > Tyr [16] 
EBC: lim Apalr=co,4) =0 [17] 


The coupling conditions between regions 1 and 2 (CCR) are given by the flow rate and pressure equalities 
at the water front radius in each layer j (LF Bittencourt Neto et al., 2020a; Nie et al., 2011): 


Api (r = rir, t) = Apyo(r = rir, t) 
Apo (r = rar,t) = Apoo(r = rar, t) 


AP (r =TnF, t) = APya(r =TnF, t) 


COR=< Xn [18] 
411 = 412 
dai = G22 
Ini In2 
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Using Darcy’s law, the flow rate CCRs can be replaced: 


CCR= 


Api (r = rir, t) = Apyo(r = rir, t) 
Apo (rT = Tap, t) = Apoo(rT = rap, t) 


AD, (7 =TlnF; t) = APya(r =TnF; t) 














OAR, _ Xo OAD 12 
(r Or ~ Rw Or 
T=T1F T=T1F 
OAP 21 — do (» PAB22 
(r Or =a OF 
T=ToF T=1oF 
OAP nt — Xo OAPy2 
(r Or ) -_ fo (r Or 
T=TnF T=!'nF 








[19] 


Between layers, there is also a coupling condition, the CCL. This condition is obtained considering that 
the pressure change in both layers is equal at the well and that the flow rate at the well is given by the sum 
of the flow rates of each layer (LF Bittencourt Neto et al., 2020a): 


CCL = 


Api, (r = est) = Apo, (r — fags) 
Apo (T = Tw, t) = ADs (7 = Tw, t) 


ADP (n—1)1(7 = Tw t) = APpi(? = Tw, t) 


Once again, the flow rate equation can be rewritten using Darcy’s law: 


CCL = 


n OAD; 





= 
Adi (Tr = Tw, t) = Ado (rT = Tw, t) 
Apo (T = Tw, t) = ADs, (rT = Tw, t) 





q = 
Qrury (AL, hi+A2,ho+...+A7, hn) > Or i 


J 


AD (n—1)1 (7 = Tw, t) = AP (7 = Tw, t) 


[20] 


There is a general pressure solution for equations [14] and [16] given in terms of the Bessel modified 
functions (Ehlig-Economides and Joseph, 1987): 


Region 1 (water): 


AD ji = Aj Ko(roj1) + By Io(roj1) 
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[22| 
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Region 2 (oil): 


AD 52 = A} Ko(ro2) + BF? Io(roj2) [23] 


For 7 =1,...,n. 
Applying the EBC in the general solution given for region 2: 


im [Aj Ko(oj2r) + B3Io(o;2")| = 0 [24] 


Knowing the properties for Bessel’s modified functions, it is true that: 





Jim | B}Jo(aj2r)] = 0 BP =0 [25] 





Hence, the pressure solution for region 2 is given by: 


ADj2 = A’ Ko(aj27) [26] 


Calculating the laplacian of both solutions: 


V?ADj1 = 05 (AjKo(ro;1) + ByIp(roj1)) = oj AD; [27] 


V?ADj2 = 05( AF Ko(raj2)) = 7j2ADjo [28] 


Replacing the equivalences above into the ODEs [14] and [16]: 


(Kj0%; — wyCepje — Xj—-1 — Xj) ADyg + Xj rAAPG—1yi — XjAVG41y: = 0 [29] 


This is a homogeneous linear system for each region i where the nontrivial solution is wanted, that is, 
Ap,; #0. That is true only if each matrix below is singular, and that implies that its determinant must 
vanish: 
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qs fork =j-—1;7>1, 


a RiP ji — WCE — Xiy = X}, for k = J, 130] 
i Xi, fork=j+1;j<n, 
0, fork Aj-—1,j7,0r 7 +1. 


For 7 = 1,2. 

Now, with the values of oj; calculated by letting det(a', ,) = 0, the coefficients of the pressure solution 
must be found in order to calculate the pressure variation. The coupling conditions are used for that end. 
Replacing equations [22] and [26] in each CCR and CCL the following equations are obtained: 

Pressure CCL: 

Aj_1Ko(oG-1)1Tw) + By_Io(oG—11Tw) = Aj Ko(oj1Tw) + ByIo(oy1Tw) [31] 
For 7 = 2,...,n. 
Flow rate CCL: 


q 





011(At Ki (rwoir) — BL (rwou)) +... + On (AL Ki (rwont) — BLT (ryon1)) = nm .. [32] 
Qty >) Ayhy 
j=l 
Pressure CCR: 
Aj Ko(oj1751) + By Ip(oq1r51) = A} Ko(oj27j1) [33] 
For j = 1,...,n. 
Flow rate CCR: 
1 1 do Ojo] y2 
A; Ki (rj ro;w) = Bil (rjrojw) = = ——|Aj Ki(rj rojo) [34] 
Aw Cjw 


For j = 1,...,n. 


Relations [31], [32], [83] and [34] provide a 3n equations that form a linear system to calculate the 
pressure coefficients Aj, Bj and A: 


q 

A} Qrurw = Muh 
By a 
a 0 
Ay 
Bi 0 
22|\=M4 0 

: 0 
Al 
Bl 0 

n 
A2 0 

i‘ 0 
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where M is formed using equations [31] to [34], and finally, it is possible to calculate the pressure variation 
at the wellbore: 


AD wf _ Al Ko(ourw) + BlIo(ourw) [35] 


Pressure profile is given by equation [35] in the Laplace field. Using the Stehfest (1970) Algorithm, this 
solution is given in the real field. 


4. Results and Discussions 


The results for the wellbore pressure profile for the analytical model proposed is presented in this 
section. A set of cases reflecting the model are presented and the accuracy of the proposed solution was 
testified by comparing it with a commercial finite difference-based flow simulator CMG (2010). In addition, 
an estimation of the equivalent permeability was presented and compared it to its real value for all cases. 

As an input for the numerical simulation, a radial grid was considered. The grid is more refined in the 
region closest to the wellbore, which is the region most affected by its presence. The oil model used in the 
simulator was blackoil. For all cases, formation crossflow is considered between the adjacent layers. 

The following parameters were considered: 


e 4 days (96 hours) injection period; 

e A piston-like flow; 

¢ Injection flow rate was defined as 500 m?/day (5.79 x 10° m?/s); 

¢ Initial pressure was defined as 300 kef/ cm?; 

e The wellbore radius was considered to be ry = 0.108 m in all cases; 

e Water viscosity will always be considered to be uy = 0.5 cP; 

e The vertical permeability was considered to be equal to the horizontal permeability for all cases; 
- SJ. and si, are considered to be 0 and Loe is calculated as in eq.(6); 


e The thickness of the shale between the layers (Ah) is considered to be zero in all cases for both the 
analytical and the numerical models. 


In order to best evaluate the efficiency of the model, the mobility ratio is considered as a parameter for 
all cases in each layer, and it is given by (LF Bittencourt Neto et al., 2020b): 


ae: 
M; == [36] 
a5 


The mobility ratio, is directly connected to the sweep efficiency of a reservoir. Low mobility ratios 
(M <1) imply in better sweep efficiency, as oil flows more easily throughout the reservoir. On the other 
hand, high mobility ratios lead to uneven fronts, reducing the swept reservoir area (LF Bittencourt Neto 
et al., 2020b). 
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Other relevant parameters can be found in table 1 for all cases considered: 









































Case | Properties | Layer 1 Layer 2 | Layer3 
k(mD) 500 1000 - 
h(m) 15 15 - 
A 
[4o(CP) 5.1 - 
M 10.2 10.2 
k(mD) 500 600 700 
a h(m) 15 15 15 
[o(cP) 0.1 
M 0.2 0.2 0.2 
k(mD) 1000 600 300 
A h(m) 15 15 15 
[o(cP) 3.5 
M 7 | 7 7 

















Table 1 — Analyzed Cases 


In addition, in the graphs of all cases, besides the pressure change curves for the analytical and numerical 
solutions, the curves of pressure derivatives with respect to the logarithm of time is also present. It is 
important to analyze the behavior of the pressure derivative as well, in order to properly interpret the 
results of the test. The curve composed by blue circles denoted the analytical solution for pressure variation 
and the curve composed by red circles denotes the analytical solution for pressure derivative. For the 
numerical solutions, yellow stars curve represents the pressure variations and purple stars the derivative 
curve. 

First, case A is analyzed. For this case, a two layer model with different values of permeability in each 
layer is considered. Pressure and pressure derivative curves for the analytical model solution and numerical 
solution are presented in figure 3 for case A: 





1 
10 Apwf (analytical) 


dApwf (analytical) 
Apwf (simulated) 
dApwf (simulated) 








** OO 








Apwf and dApwf (kgf/cm2) 
3 














Fig. 3 — Analytical and numerical pressure variation and derivative solutions for case A 


Even though there are no different regions of permeability in the layers for this case, the derivative curves 
behaves as if there were. That is because, the regions of equal permeability are filled with different fluids. 
The oil and water viscosity values for this case, causes the value of the mobility ratio M to be grater than 1, 
implying on higher pressure derivative values on early times. It is possible to see a close agreement between 
the analytical and numerical simulated curves. 
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Case B considers a value of oil viscosity of wu, = 0.1 cP, which was 5 times smaller than the water’s. 
Implying on a value lower than 1 for M, in this case it was 0.2. This case also considers different values of 


permeability per layer. The graphs for pressure and pressure derivative for case B are represented in figure 
4: 





Apwf (analytical) 
dApwt (analytical) 
Apwf (simulated) 
dApwt (simulated) 


** OO 











Apwf and dApwf (kgf/cm2) 














10°? : : 
10° 10° 10° 
Time (h) 


Fig. 4 — Analytical and numerical pressure variation and derivative solutions for case B 


It is possible to see that the values in the derivative curves at initial times were five times smaller than in 
the final times. The derivative curves were increasing because M <1and the sweep efficiency was improved. 
There was a close agreement between the numerical and analytical curves for this case also. 

Finally, case C' is considered. That case considers three layers with different properties and fig = 3.5 cP. 
The graphs for pressure and pressure derivative for this case are presented in figure 5: 





1 
10 Apwf (analytical) 


dApwf (analytical) 
Apwf (simulated) 
dApwt (simulated) 


** OO 











Apwf and dApwf (kgf/cm2) 














Fig. 5 — Analytical and numerical pressure variation and derivative solutions for case C’ 


Similarly to case A, case C’ has M> 1, however, since the difference between the fluids viscosity values 
is not as big as in case A, the difference in the values of the derivative curves is not as great. There was 
a small, and apparently constant, difference between the analytical and numerical pressure curves. One 
reason for that might be because there is a big difference in the permeabilities values from one layer to 
another causing a bigger awareness of the presence of formation crossflow in the reservoir, and the way the 
simulator includes the crossflow is different than the way presented in this work. However, the behavior of 
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the curves is still very similar, resulting in a fair approximation. 


A. Equivalent Permeability. The pressure measured at the well is directly linked to the equivalent perme- 
ability (Keq). In that sense, another way to evaluating the efficiency of the model presented in this work is 
by estimating the equivalent permeability. 

For a multilayered reservoir, like presented in Cobb et al. (1972), it is given by: 


[37] 


Since the model presented here reduces the cases of different regions radii into one with regions of equal 
radii, then the derivative profile enables the computation of equivalent permeabilities that combine the 
respective regions of the layers. That is, to obtain the first equivalent permeability, region 1 of layer 1 is 
combined with region 1 of layer 2 and so on. 

Using the source line logarithmic approximation, the following equation is used to estimate the equivalent 
permeability, derived from the single phase formulation(Cobb et al., 1972): 


>= _ Ho 
“1 Qhrm 2 


= 





where m is the constant derivative level and hr is the total thickness. 
In table 2, the estimated and real equivalent permeabilities for the region filled with water, are presented 
for the cases considered in this chapter along with the percentage errors of the approximation: 





keq estimated (mD) | Error (%) 
748.24 0.49 
600.33 0.01 
629.64 1.48 

















Table 2 — Real and estimated equivalent permeability values and percentage error for cases A to C 


The cases presented errors less than 2%. This indicates that the proposed formulation may be useful in 
obtaining reservoir parameters. 


5. Conclusions and Suggestions for Future Work 


To the best of the authors” knowledge, this model to solve the injectivity test pressure response problem 
using a radial composite approach, considering formation crossflow, has never been presented in literature 
before. 

Based on the work presented by Ehlig-Economides and Joseph (1987) which considered formation 
crossflow, the present work extended it for a two phase model using a composed reservoir approach. 

For all cases, a close agreement between the approximated analytical solution and numerical data 
was achieved. Results show that the suggested analytical model is able to describe pressure behavior in 
multilayered reservoirs during injectivity tests using a radial composite approach. Furthermore, results 
show a small error when comparing the estimated equivalent permeability to the real ones for all cases 
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presented. This work extends the comprehension of pressure behavior to a wider range of scenarios and the 
results obtained might be used to reservoir management applications. 

The presented solution considers a zero skin factor. A solution which considers different values of the skin 
factor is of practical interest. In addition, this work considers a single region of permeability in each layer. 
A solution which considers multiple regions of permeability per layer is also of practical interest. For the 
implementation of that last suggestion, there would be more semipermeability coefficients per adjacency. 


6. Nomenclature 


ajz - Matrix elements 

Aj, Bj - Coefficients for jth layer and ith region defined in the pressure solution 
ce - Total system compressibility 

hj - Thickness of layer j 

hr - Total reservoir thickness 

I;, K;, - Modified Bessel functions of first and second kind and order 7 

Keg - Reservoir equivalent permeability 

ki - Horizontal permeability in layer j and region i 

n - Number of layers in reservoir system 

m, - Number of regions in layer j 

p; - Initial Pressure 

pi - Reservoir pressure in layer j and region i 

Ap - Pressure change pj - Pressure change in the Laplace domain 

Pw - Well-bottom hole pressure 

qji - Flow rate for layer j and region i 

r - Radius 

ri - Radius of region i in layer j 

re - Reservoir outer radius 

Tw - Wellbore radius 

t - Time 

u - Laplace variable 

Xi - Semipermeability coefficient between layers 7 and j + 1 in region i 
xi - Vertical permeability-thickness ratio 

t - Permeability-thickness product in layer j and region i 
@; - Porosity in layer j 

[up - Fluid viscosity 

j1;; - Viscosity of region i in in layer j 

w; - Porosity-thickness product in layer j 


K 
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